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Planetary rings are common in the outer Solar System but their origin and long-term evolution is still a matter of 
debate. It is well known that viscous spreading is a major evolutionary process for rings, as it globally redistributes the 
disk's mass and angular momentum, and can lead to the disk's loosing mass by infall onto the planet or through the 
Roche limit. However, describing this process is highly dependent on the model used for the viscosity. In this paper 
we investigate the global and long-term viscous evolution of a circumplanetary disk. We have d eveloped a simple ID 
numerical code, but we use a physically realistic viscosity model derived from N-body simulations ( Daisaka et all l200lh . 
and dependent on the disk's local properties (surface mass density, particle size, distance to the planet). Particularly, 
we include the effects of gravitational instabilities (wakes) that importantly enhance the disk's viscosity. This method 
allows to study the global evolution of the disk over the age of the Solar System. 

Common estimates of the disk's spreading time-scales with constant viscosity significantly underestimate the rings' 
lifetime. We show that, with a realistic viscosity model, an initially narrow ring undergoes two successive evolutionary 
stages: (1) a transient rapid spreading when the disk is self-gravitating, with the formation of a density peak inward 
and an outer region marginally gravitationally stable, and with an emptying time-scale proportional to 1/Mq (where 
Mo is the disk's initial mass) (2) an asymptotic regime where the spreading rate continuously slows down as larger parts 
of the disk become not-self-gravitating due to the decrease of the surface density, until the disk becomes completely 
not-self-gravitating. At this point its evolution dramatically slows down, with an emptying time-scale proportional to 
1/Mq, which significantly increases the disk's lifetime compared to the case with constant viscosity. We show also that 
the disk's width scales like t 1 / 4 with the realistic viscosity model, while it scales like t 1 ^ 2 in the case of constant viscosity, 
resulting in much larger evolutionary time-scales in our model. We find however that the present shape of Saturn's rings 
looks like a 100 million-years old disk in our simulations. Concerning Jupiter's, Uranus' and Neptune's rings that are 
faint today, it is not likely that they were much more massive in the past and lost most of their mass due to viscous 
spreading alone. 
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1. Introduction 

Saturn's rings are one of the most puzzling object of 
the Solar System. They've been studied for over 400 years 
since their discovery by Galileo in the early 17th century, 
but still many key questions remain unanswered. In par- 
ticular, how old the rings are is still a matter of debate. 
While published scenarios for their origin suggest an early 
formation in the history of the Solar System : (1) rem- 
nant of Saturn's sub-nebula disk dPollack et al" . 19731). (2) 



tidal d isruption of a comet (jDonesl . 1 19911 : ICharnoz et al 



2009 3), (3) destruction of a satellit e insid e Saturn's Roche 
zone (Pollack et al.. 1973: Harris . 1984: Charnoz et al 



200911) . observations and theoretical arguments such as 
the viscou s spreading of t he A ring in a few hundred mil- 
lion years (lEspositol.ll986f). or the low m eteoritic pollution 
of the rings (jCuzzi and Estrada! . I1998T ). lead to the con- 
clusion that they must be quite young. Thus we are in a 
paradoxical situation. 

The rings' evolution is regulated through 3 main phys- 
ical processes. (1) Viscous spreading: particle collisions 
dissipate energy while conserving the total angular mo- 
mentum, resulting in the spreading of the disk with the 
mass bein g transferred inward and t h e ang ular momentum 
outward dLvnden-Bell and Pringlel. 19741: Goldreich and 
Trcmaine, ll982T ). (2 ) Interactions with the planet's satel- 
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lites at resonances (| Goldreich and Trcmain q. | l979f). an d 
(3) Meteoritic bombardment (jCuzzi and Estrada 



1998). 



All these processes modify the disk's angular momentum 
and transport material, which can thus leave the disk when 
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it reaches the planet's radius or crosses the Roche limit 
where it can start to be accreted into small moons (Charnoz 
et al., In press). These processes can then cause a drop 
of the disk's mass. Seeing how they strongly influence the 
possible lifetime of the rings, the paradox stated above 
must come from an insufficient understanding of these pro- 
cesses on the long term and from a lack of models. A more 
detailed and global study of these processes over long time- 
scales is needed. 

While process (3) is specific to Saturn's rings, processes 
(1) and (2) are common to other disks in the universe, for 
instance protoplanetary disks, accretion disks, or galax- 
ies. But Saturn's rings are particular in the sense that 
they are devoid of gas and composed only of macroscopic 
particles, so that pressure and radiative effects can be ne- 
glected compared to gravitational effects. Consequently, 
while the approach remains valid, results on the evolution 
of a gaseous disk under the aforementioned processes can- 
not be directly transposed to Saturn's rings but have to 
be studied considering their specificities. 

The viscous spreading of a di sk has already been stud- 
ied, particularly accretion disks j Lvnden- Bell and Pringld . 
1974 : iBath and Pringld . 1 198 it IPringld . Il98lh . However, 
this process is highly dependent on the model used to de- 
scribe the disk's viscosity. Several important works have 
been publish ed on local properties of the rings ( Wisdom 
and Tre maine ll988l:ISald.ll992l:lRichardson).ll994l: Daisaka , 
and I da, 1 19991 : lOhtsuki and Emoril . l2000l : iDaisaka et all 
20011 ) , resulting in precise estimates of the viscosity depen- 
dence on the disk's parameters (surface density, particle 
properties, ...), but these local properties have not been 
studied yet on large spatial and temporal scales. In this 
paper we focus our study on the viscous spreading of the 
disk considering a physically realistic viscosity model and 
the evolution of the whole disk. Eventually, constrain the 
evolution of the rings will require to consider many dif- 
ferent physical processes, but this is for now well beyond 
current computer capacity. Moreover, questions such as 
the spreading rate of the rings, the amount of material 
falling onto the planet or leaving the rings through their 
outer limit, are important issues that still need to be ad- 
dressed. 

Published estimates of the rings' viscous age are vis- 
cous time-scales Ar 2 /i/ that give the time needed for a 
ring with constant viscosity v to reach the width Ar. How- 
ever, we can expect the viscosity of the rings to be a rapidly 
varying function of the disk's local properties (surface mass 
density, ring particle's radius, distance to the planet), so 
that such time-scales might be substantially inaccurate. 

In this work we study the effects of a time-dependent 
and space-dependent viscosity model on the long term vis- 
cous evolution of the complete ring system. We perform 
1-dimensional numerical simulations, along the radial di- 
rection, of the disk's viscous evolution. This implies sim- 
plifications of the physics but allows for simulations over 
many dynamical times (typically the age of the Solar Sys- 
tem). As a consequence, our results will be a first-order 



study of the global viscous evolution of dense planetary 
rings. However, it may give us clues to understand what 
has led to the rings we can observe today. Moreover, this 
approach has already been used with i mportant results 



in the study of protoplan etary disks (e.g. iTakeuchi et al 
19961 : lAlibert et all 1200.^ . 



In section 2 we describe the viscosity model and the 
numerical procedure. In section 3 we study the viscous 
evolution of an initially narrow ring over 5 billion years us- 
ing two models: constant viscosity (hereafter CV) and non 
constant viscosity (hereafter NCV), and develop analyti- 
cal time-scales to describe the disk's evolution. In section 

4 we analyse the influence of ring parameters. In section 

5 we summarize the results and discuss the evolution of 
today's rings of Saturn in the light of these results. 

2. Numerical method and viscosity model 

In this section we present the hydrodynamical equa- 
tions describing the evolution of the rings, the viscosity 
model we use, and the principles of the numerical code we 
have developed. While adopting a very general approach, 
we present our results for the case of Saturn's rings. 

2.1. Basic equations 

We use a 1-dimensional approach to study the viscous 
evolution of a Keplerian pressureless disk. To compute the 
evolution of the surface density S (R, t) of an elementary 
annulus of the rings located at distance R f rom the plane t 
and at time t, we use the same approach as IPringld fll98lh . 
The equation of mass conservation reads: 



R 
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dt + dR 



{RZv K ) = 0, 



(1) 



where vr (R, t) is the radial velocity. The angular momen- 
tum conservation reads: 



RdR\ dR dt K ' 



RdR 
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where — y/ GM/R 3 is the orbital frequency, with G the 
gravitational constant and M the planet's mass, and v is 
the kinematic viscosity. Combining Eqs. ([1]) and ([2]) and 
replacing f2 by its expression we get a single equation for 
the temporal evolution of the surface density: 



<9S 
dt 



3_d_ 
RdR 



(3) 



2.2. Viscosity model 

2.2.1. Description of the model 

The viscosity of a disk is related to angular momen- 
tum transport through particle interactions. The angu lar 
momentum flux reads ( Lvnden-Bell and Pringld . H97l l 



$ = 37r£fti?V 



(4) 
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Using the Boltzman n equation, one can der ive this angular 
momentum flux as (jTakeda and Idal . 120011 ) 



d_ 
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where Ug is the mean tangential velocity averaged az- 
imuthally. The viscosity can then be computed by esti- 
mating the angular momentum flux and inverting Eq. (j4|). 

Angular momentum transport can be divided in three 
components, each one related to a specific viscosity. The 
translational viscosity ftrans is related to the transport of 
angular momentum due to the random motion of particles 
(usually referred to as the "local" co mpone nt, see Goldre- 
ich and Tremaine (|l978l )). Due to their finite size, angular 
momentum is also transported via sound waves travelling 
between the centres of colliding particles (usually referred 
to as the "non-local" compon ent). This is represen t ed by 
the collisional viscosity p cn u (Araki and Tremaine] . 1986; 
Wisdom and Tremainel . Il988f ). Finally, angular momen- 



tum is transported by gravitational scatteri ng of parti- 
cles due to the presence of self-gravity wa kes jgalol . (l992l . 
19951 : iRichardsonl. Il994l: iDaisaka and Idal. Il999l: Ohtsuki 



and Emori. l2000f ) . This is represented by the gravitational 
viscosity ^g rav . 

The self-gravity wakes are gravitational aggregates in- 
duced by the effects of both self-gravity and collisional 
damping of particles. Outside of a wake, particles move 
randomly, but inside their motion becomes coherent, which 
yield s systematic motion w ith large bulk viscosity (jSalol . 
ll995l : lDaisaka and Idal . [l999] ). The wakes modify the angu- 
lar momentum transport , and thus the viscosity, through 
the g ravitational torqu e ( Larsonl . I1984T ) and the wake mo- 
tion ( Lin and Pringlel . 1987t ). Two regimes must then be 
considered, whether the disk is gravitationally stable (no 
wake s) or not . This is measured by the Toomre Q param- 
eter (jToomrd . Il964r ) 



Q = 



3.36GS' 



where cr r is the particle radial velocity dispersion. Even 
though Toomre showed that the disk is gr avitat i onally 



though loomre snowed that the disk is gr avitat i onally 
unstable for Q < 1, N-body simulations by ISalol ( 1995h 
showed that wakes start to form for Q < 2. 

In our simulations, we use the following formulations 
for the viscosity components. We set the transition be- 
tween the two regimes, self-gravitating and not-self-gravi- 
tating (hereafter SG and NSG), strictly at Q = 2. 



^coll 



x grav > 



with 



^trans 
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[26r 



if Q > 2, 
if Q < 2, 



; 2 s 2 



fcoii = r p flT V Q, 







if Q > 2, 



^trans if Q < 2, 
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(8) 

(9) 
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where r = 3S/ (4r p/ o p ) is the optical depth with r p and 
p p the particle's radius and density, is a dimension- 
less parameter equal to the ratio of the particle's Hill 



(5) radius rh to its physical diameter: 



rh/(2r p ) with 



1 /3 

rh = (2TO p /(3Af;s)) ' R, and m p is the particle's mass. 
It's basically a dimensionless measure of the distance to 
the planet, while the optical depth serves as a dimension- 
less measure of the disk's surface mass density. 

In the NSG regime (Q > 2), the translational vi s cosity 
is an analytical result from lGoldreich and Tremainel ( 1978h 
(Eq. dSJ, top) , and the collisional v i scosity is an analyti- 
cal result from lAraki and Tremaine ( 1986 ) (Eq. ^ ) . In 
the SG regime (Q < 2) (Eqs. [51 bot t om fc UM bottom), 



we use the results of Daisaka et al. (2001). as they are 



the first to include the effects of the wakes in the calcula- 
tion of the viscosity. They performed multiple 3D shear- 
ing box N-body simulations, including the effects of self- 
gravity wakes, with lm-radius particles. Their set of pa- 
rameters includes high and low densities, which is suitable 
to track the rings' evolution in a large variety of physical 
conditions. Note that although the collisional viscosity is 
also enhanced in the presence of wakes, its value is about 
one order of magnitude smaller than the translational and 
gravitational visc osity in that case, ac cording to the sim- 



ulation results of IDaisaka et al.l (|200lf ). Thus we use the 



same prescription for the collisional viscosity whether the 
disk is gravitationally stable or not. 

Th e parameter is also expressed in IDaisaka et al 
(2Q02) as |F grav |/|F coll | = 4r* 2 , where F grav is the self- 



gravity force between particles, and F co \\ is the impulsive 
force exerted on particles during collisions. This can be 
used to evaluate the velocity dispersion er r . If rjj < 0.5, 
-Fcoii dominates and the velocity dispersion is the relative 
Keplerian velocity between particles <r r = 2r p f2. Con- 
versely, if r£ > 0.5, -Fgrav dominates and the velocity dis- 

:a pe ve- 
al. Il999l : 
> 0.5 so 



persion is regulated to be the particle's surface esca 
y/Gmp/rp (jSalol . Il995c IDaisaka and Id 



locity o 



(6) lOhtsukil . Il999l ). For Saturn's rings we have 



in the following we always consider that cr r = y/ Gm p /r p . 
Note that we do not consider here the effects of thermal 
conduction that could modify the velocity dispersion. 

Toomre's Q parameter depends on the disk's local prop- 
erties: it decreases with surface mass density E and dis- 
tance to the planet R. Some parts of the disk may then 
be gravitationally stable while others are not, and we can 
expect theses regions to undergo distinct evolutions. As 
the disk evolves, the surface density will decrease because 
of the disk's spreading, and the Toomre Q parameter will 
increase. Then regions initially self-gravitating may be- 
come not-self-gravitating as the disk evolves, leading to 
local modifications of the viscosity. Note that with the 
adopted prescription of the velocity dispersion, Q is also 
proportional to r p . So a disk with bigger particles will be 
more likely to become not-self-gravitating. 
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2.2.2. Viscosity dependence on ring parameters 

The viscosity components depend on several parame- 
ters that are likely to vary within the rings: particle radius, 
surface mass density and distance to the planet (Eqs. flSJ) 
and ©). Even though this model may be used for any 
planetary ring, we are most interested in Saturn's rings. 
Here we discuss the resulting viscosity for the specific case 
of Saturn. As an example, we plot in Fig. [T]the variation 
of t'trans and i^coii, f° r a constant value of the surface den- 
sity (400 kg.m~ 2 , the surface density of today's A ring, 
Tiscareno et all 120071 ). with respect to R and for different 



values of r p . 
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Fig. 1. Radial evolution of the viscosity for different ring par- 
ticle sizes. The surface density is set to 400 kg.m -2 . The trans- 
lational viscosity (black curves) increases with the distance to 
the planet, while the collisional viscosity (red curves) decreases 
with the distance. The translational viscosity can be several 
orders of magnitude larger in the self-gravitating case (black 
dotted line). 

Both translational and collisional viscosities increase 
with particle radius, but they have opposite dependences 
on the distance to the planet: the collisional (translational) 
viscosity decreases (increases) with R. In the SG regime, 
the particle radius does not modify the value of the vis- 
cosity itself, but it impacts Toomre's Q parameter which 
controls the areas of the rings that are in the SG regime. 

As the disks spreads, its surface density will decrease. 
Transitions from SG to NSG regimes might thus occur, 
leading to important variations of the viscosity. To iden- 
tify the effect of a self-gravity regime transition on the 
viscosity, we study the relative magnitude of the theoreti- 
cal viscosities in the two regimes, and analyse the position 
of the transition Q = 2 in the space of parameters (R, £). 
The result is plotted in Fig. [2j for particles with r p = I m. 

A transition between the SG to the NSG regime will 
occur at a given distance R when the surface density £ 
goes below the limit represented by the dashed-black curve 
for which Q = 2. Close to the planet, the transition from 
the pink to the dark blue region will lead to an increase 
of the viscosity, because in these regions vsg < ^nsg- Far 




100 120 
Distance to Saturn (1000 km) 



Fig. 2. Relative magnitude of the total viscosity i^trans + 
^coii + fgrav m the self-gravitating and not-self-gravitating 
regimes, and position of the regime transition. The dashed 
curves are the transition between the self-gravitating and not- 
self-gravitating regimes at Q = 2, and the set of parameters 
(R, £) for which the theoretical viscosities in the two regimes 
are equal. 



from the planet, the transition from the green to the yellow 
region will lead to a decrease of the viscosity, because in 
these regions z^sg > ^nsg- The transitions between the 
two behaviours occurs at R ~ 78, 000 km. 

This might be explained by looking into the particle 
velocity dispersion. Inside a wake, this velocity is the es- 
cape velocity: <r r = */ Gm p /r p . In the absence of wakes, 
the velocity dispersion is the difference of orbital speed be- 
tween two particle distant by 2r p : cr r = 2r p i7. Far from 
the planet, the velocity in a wake is greater, while close 
to the planet the difference of orbital velocity dominates. 
Equalizing the two expressions shows a transition around 
R r-j 84, 000 km, which is in good agreement with the tran- 
sition observed in Fig. 

2.3. A simple code 

2.3.1. Numerical procedure 

We have developed a ID finite elements code on a stag- 
gered mesh: the surface density is estimated on a regular 
grid, and its evolution is computed by estimating the mass 
fluxes at the inner and outer edges of each bin of the grid. 
This allows a second order accuracy in space derivatives. 

By integrating Eq. ^ over the bin width dR, the 
variation over a time step dt of the surface density in the 
bin i is given by: 



Flux (Ri + ^) - Flux (Ri-tf) 
Bin Surface Area 



The flux at time t reads: 



Flux(i?) = 6tt 



R M (^) 



dt, (ii) 



(12) 
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It is estimated using a second-order (in space and time) 
explicit Runge-Kutta scheme. 

As in any numerical procedure we can't study scales 
smaller than the grid resolution (~ 100 km per bin). As a 
result we don't expect to see small-scale instabilities such 
as self- gravity wakes, or viscous overstabilities ( Latter and 
Ogilvie, J00.9J) , appear in our simulations. Moreover, as we 
use a constant value for the velocity dispersion, we do not 
satisfy the conditions for viscous instabilities to appear. 
Anyway the purpose of this work is to study the rings evo- 
lution on large scales, and even though the full resolution 
of the hydrodyna mics equations in 2 or 3D wo uld be desir- 
able ( e.g. like in Baruteau and Massed. 2008 ; Masset and 
Casoli, |2009j)~ such codes arc limited to short time-scales. 
We keep the main physical ingredients in this ID model, 
allowing us to perform simulations over 5 billion years. 

2.3.2. Validation 

To validate numerically the behaviour of our code we 
check that the system's mass is numerically conserved. We 
perform a simulation starting with an initially narrow ring 
with Gaussian shape and an arbitrary initial mass. The 
boundary conditions are set so that material cannot pass 
through the limits of the domain: the disk is "trapped" 
inside the domain. We let the ring evolve viscously and 
compute the total mass of the disk using the output surface 
density. The initial viscous time-scale AR 2 /v is ~ 11, 500 
years for AR = 3000 km. After 100 million years of evolu- 
tion (about 9000 viscous time-scales), the mass is equal to 
its original value down to the machine precision (10~ 18 ). 
Our code numerically conserves the system's mass. 

Equation [3] includes the equation of angular momen- 
tum conservation, so by construction our code should also 
conserves angular momentum, as long as the disk has not 
reached the domain's limits where the boundary conditions 
cause modification of the disk's angular momentum. Dur- 
ing the period where the disk spreads freely between the 
domain's limits, we have checked that angular momentum 
is also numerically conserved down to machine precision. 

In order to validate the general beh aviour o f the c ode, 
we investigated the solution given by IPringle ( 1981 ) for 
the viscous spreading of an accretion disk with constant 
viscosity. Our numerical results are in good agreement 
with the analytical solution, even though slightly different 
because the width of the radial grid does not allow to define 
an initial ring with a perfect Dirac distribution. 

3. Viscous spreading with variable viscosity 

Here we study the viscous evolution of a "standard 
model" using CV and NCV models. We present first the 
disk initial conditions, then the viscous evolution over 5 
billions years, and finally we focus on the evolution of the 
disk's mass. 



3.1. Description of the standard model 
3.1.1. Initial profile 

The initial state of Saturn's rings is dimly constrained, 
as it depends on the formation scenario, which is still sub- 
ject to discussions. For instance the scenario based on the 
destruction of a satellite via a comctary impact implies 
that the satellite that created the rings must be outside 
the synchronous orbit an d inside the Roche lim it for ice 
at the time of the impact ( Charnoz et al. . 2009bf ). But we 



have no further constraints on the initial rings width, posi- 
tion, surface mass density, particle size distribution. . . For 
the sake of simplicity and generality, the initial conditions 
are defined as follows: 

Initial surface density profile. We consider at t — that 
the disk is a narrow ring ~ 3000 km wide with Gaussian 
profile. This is useful to study its evolution over several 
viscous time-scales. The ring centre is set slightly outside 
the synchronous orbit (~ 110,000 km), which is more or 
less the middle of today's ring system. The influence of 
the initial profile is studied in section 4. 

Particle radius. The particle radius of today's rings ranges 
from a few centimetres to several metres. There are also 
radial variations, and o n the same orbit d i fferen t sizes 
of p articles may coexist dCuzzi and Pollack!. 19781: Cuzzi 



et al.. Il980l : IPorco et alll2008fi. Particle size s shou ld also 
evolve through collisions (jAlbers and Spahn . 2006). Un 



fortu nately the viscosity model developed by Da isaka et al 



(2001) does not include a particle size distri bution, but a 
fixed particle radius r p . iMarouf et al. I dl983h showed that 
particles in Saturn's rings follow a power law iV(r D ) oc r~' 



Shu and Steward (|l985h 



with a sa 3 for R\ < r p < i?2- 
then computed that this distribution can be represented 
by a single equivalent particle radius R e given by R e = 
\/3Ri/it. For Ro ~ 10 m we get R„ ~ 5 m. Goldreich and 
Trcmaine (Il982h suggests an equivalent radius r p = 1 m, 
which is also a common value in many N-body simula- 
tions. Considering the uncertainty around an equivalent 
particle size, we choose r p = 1 m for our standard model, 
and discuss the effect of particle size in section 4. 

Initial disk mass. Considering that today's rings mass i s 
estimated to be about Mimas' mass (jEsposito et al.l . ll983T ). 
one can assume that the initial disk mass was at least Mi- 
mas' mass. But on the other hand meteoritic bombard- 
ment is believed to bring a larg e amount of mass into the 
rings (jCuzzi and Estra da. 1998), which could substantially 



modify the rings' mass. For these reasons we use in our 
standard model an initial mass equal to Mimas' mass, and 
study the initial mass influence in section 3.3. 

3.1.2. Constant viscosity values 

To analyse our results with the NCV model, we also 
perform simulations with a constant and uniform viscosity 
model. Assuming that the rings were denser in the past, 
which would increase their viscosity, we choose a value 
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Fig. 3. Disk surface density at different evolution times for variable (solid line) and constant (dashed line) viscosities, (a) Initial 
profile, (b) At 1000 years of evolution, (c) At 10 4 years of evolution, (d) At 10 5 years of evolution. The disk with variable 
viscosity spreads faster and does not keep the original shape of the density profile. 



of 0.1 m 2 .s 1 , which is about ten tim es the viscosity of 



today's A ring (jTiscareno et all [2007). Anyway, in the 



CV case, the value of the viscosity does not modify the 
density profile, but only the viscous time-scales. 

3.1.3. Boundary conditions 

The radial domain ranges from R — 65, 000 km to 
R = 140,000 km (R = is the centre of Saturn). The 
boundary conditions are set as follows: material is released 
from the disk when it passes the domain's limits. For the 
inner boundary we consider that the material falls onto 
the planet and is destroyed, and for the outer boundary 
we consider that the material crosses the Roche limit and 
starts to be accreted in small moons that do not inter- 
act viscously with the disk. They should however interact 
with the disk at resonances, but this is be yond the scope 
of this paper. This has been investigated bv lCharnoz et al. 
(In press). 

We remove the material passing through the domain's 
limits, and the corresponding angular momentum, by set- 
ting the surface density of the first and last bins to zero at 
t = 0, and at each time step we equalize the inward and 
outward fluxes in these bins. As a result, the total flux on 
the first and last bins is 0, and all material entering these 
bins is immediately removed, along with the correspond- 
ing angular momentum. No torque is exerted on the ring 
by the boundary. 

3.2. Viscous spreading over 5 billion years 

In this section we study the evolution of the standard 
model, using constant and variable viscosities, over the age 
of the Solar System (~ 5 billion years). 



3.2.1. Relevant parameters 

To analyse our results we study the evolution of three 
quantities: 

• the surface density profile; 

• the total mass of the disk, derived from the surface 
density, in order to track the mass lost during the 
viscous spreading; 

• the total mass passing through the left and right lim- 
its of the domain. 

3.2.2. Early evolution (0 - 10 5 years) 

In Fig. [3] the surface density at different times for both 
CV and NCV models is plotted: at and after 10 3 , 10 4 
and 10 5 years of evolution. As in section 2.3.2, the initial 
spreading time-scale is ~ 11, 500 years. 

In the CV case, the disk edges remain smooth during 
the spreading, while steep edges are created in the NCV 
case (Fig. \3]p). In this latter case, v being an increasing 
function of S, the viscosity is high in the centre and low at 
the edges. Thus the centre spreads faster than the edges. 
As a result, a lot of material is transported from the core 
towards the edges, but the edge itself fails to spread this 
great amount of material and is "overwhelmed" by mate- 
rial coming from the core (Fig. [3Jd, solid line). As a result, 
the surface density rapidly increases close to the edges, 
which become steep. 

At the beginning the NCV disk appears to spread faster 
than the CV disk: after 10 5 years of evolution (Fig. |3Ji), 
the disk width is about 1.5 x 10 4 km in the NCV case, 
and 10 4 km in the CV case. This can be explained by the 
fact that the initial value for the CV model is too small. 
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Fig. 4. Disk surface density at different evolution times with variable (solid line) and constant (dashed line) viscosities, (a) At 
lMyrs of evolution, (b) At f 00 Myrs of evolution, (c) At 1 Gyrs of evolution, (d) At 5 Gyrs of evolution. The disk with constant 
viscosity is emptied in ~ 10 9 years, while the disk with variable viscosity remains massive over 5 Gyrs with a density peak inward 
and lower densities outward. 



The computed viscosity at the rings centre is initially ~ 
100 m 2 .s _1 , much larger than the 0.1 m 2 .s _1 value we 
use for the CV model. But after 10 5 years of evolution 
the viscosity drops to ~ 0.8 m 2 .s -1 . Thus whereas the 
initial spreading of the NCV disk is very rapid due to high 
viscosity, as soon as the disk has spread, the viscosity drops 
and the evolution slows down dramatically. 

3.2.3. Evolution over 5 billion years 

In Fig. @]the surface mass density at 10 6 , 10 8 , 10 9 and 
5 x 10 9 years is plotted. 

We see here that the CV disk has spread much faster 
than the NCV disk (Fig. 2k). After 1 billion years, the 
CV disk is no longer visible because almost all material 
has been spread out of the rings (Fig. 0] (c) and (d)). 
Its average surface density is 5 x 10~ 5 kg.m~ 2 . The NCV 
disk has evolved significantly, and memory of the Gaussian 
initial conditions has been completely lost. 

The NCV disk evolves only very little from 1 to 5 billion 
years. This is due to a drop of the viscosity because of the 
surface mass density being now much lower than initially, 
resulting in a much slower evolution. After 1 billion years 
of evolution, the rings viscosity is ~ 10~ 3 m 2 .s _1 . 

In conclusion, this study illustrates that spreading time- 
scal es estimated using only A R 2 / v with a constant y iscos- 
ity ( Goldreich and Tremainel . 1982 ; Esposito . 198(f ). are 
quite inaccurate when compared with our numerical simu- 
lations that include a realistic viscosity model. As a result 
the viscous spreading rate of the rings is far from being 
linear. 



3.2.4- Evolution of the disk's edges 

Initially, due to the Gaussian shape we use, the disk's 
edges are smooth. Due to a low surface density, the region 
close to the edges is initially in the non-self-gravitating 
regime. When the disk starts to spread, material is trans- 
ported from the core toward the less viscous edges. This 
increases the surface density close to the edges, which then 
become self-gravitating. As the disk continues to spread, 
the surface density close to the edges decreases. Eventu- 
ally, the edges become non-self-gravitating again. 

As mentioned in section 2, the transitions from SG to 
NSG areas should lead to changes in the spreading of the 
disk, with different behaviours depending on the distance 
to the planet. Close to the planet, this leads to an increase 
of the viscosity (see Fig. Q] black dashed and dotted lines), 
which becomes inversely proportional to R: the closer to 
the planet, the higher the viscosity (as it is in this case 
dominated by the collisional component). As a result, the 
inner edge is smoothed (Fig. [5j top). 

Far from the planet, the opposite process occurs. Here, 
the transition from SG to NSG results in the viscosity 
dropping by a factor ^100. The material flux is propor- 
tional to V(i/E) (Eq. ([21)). When the viscosity is constant, 
this can be rewritten as i/VS, so that the flux depends only 
on the surface density gradient. With variable viscosity, 
this simplification is no longer valid: the flux depends on 
the gradient of both the viscosity and the surface density: 
V(z/E) = (Vi/)E + i/(VE). 

Firstly, in the NSG regime, the viscosity dependence 
on R is t'trans oc R?l 2 and v co \\ cx i? -3 / 2 , while in the SG 
regime v w J-'trans oc i? 19 / 2 (Eqs [5] and |5J) . As a result, 
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Fig. 5. Zoom on the disk's edges close to the domain's limits. 
The inner edge is smooth and progressive (Top), while the outer 
edge is sharp (bottom). 

the viscosity gradient is much lower in the NSG regime 
than in the SG regime. Secondly, due to the presence of a 
steep edge, the density gradient is larger close to the edge 
than in the core. However, when the edge enters the NSG 
regime, the viscosity significantly decreases and damps the 
large surface density gradient. Globally, the material flux 
is greater close to the outer edge than at the edge, which 
causes material to accumulates and creates the observed 
plateau (Fig. [5j bottom). 

These opposite behaviours can be more appreciated 
using the self-gravity-regime-transition diagram presented 
in section 2.2.2. In Fig. [5] the disk's surface density at 
t = 1 Gyrs is plotted over the map of the different viscos- 
ity regimes as a function of R and E. It shows in particu- 
lar that the outer plateau closely follows the limit Q = 2. 
Note that this marginally gravitationally stable plateau is 
similar to Saturn's A ring, which is known to be slightly 
gravitationally unstable, with the re gular appearance an d 
disappearance of self-gravity wakes (jFerrari et all l2009f ) . 

3.2.5. Spreading time-scales 

We look for semi-analytical solutions of Eq. (J3j) to ad- 
dress the difference of spreading time-scales between the 
disk with constant and variable viscosities. In a classi- 
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Fig. 6. Disk's surface density at t = 5 Gyrs (solid black 
curve), and self-gravitating regime. Regions of the disk for 
which the surface density is above the dashed-curve "Q = 2" 
are self-gravitating, while those below are not-self-gravitating. 

cal mass diffusion process the standard deviation of the 
mass distribution AR 2 , increases linearly with time, so 
that AR 2 jv is proportional to time. In Fig. [7] we plot 
this quantity for the CV and NCV disks. In order to 
avoid any boundary effect, we restrict the analysis to the 
t < 5 x I0 6 years where the disk has not yet reached the 
domain's limits. 




Fig. 7. Evolution of the ratio disk's width over viscosity, 
against time, for a disk with constant (dashed line) and vari- 
able (solid line) viscosity. The diamonds and triangles are data 
points from simulations. Note how AR 2 /is remains propor- 
tional to time. 

The width of the ring is defined as 

AR 2 = E ' (^-flc) 2 ^ (13) 
where R c is the "centre" of the disk defined by 

R c = ^^i, (14) 
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For the disk with variable viscosity, we have used the mean 
viscosity, weighted by the bin's mass 



< v >= 



(15) 



where Mj = 27ri?i(ii?E J ; is the bin's mass and dR is the 
bin's width. 

It is remarkable that despite the change of viscosity 
by several orders of magnitude during the evolution of the 
system, AR 2 /v remains closely proportional to time at 
any instant (Fig. [7])- However, this expression can no 
longer be used as a "rule of thumb" to estimate, with the 
initial disk's parameters, the spreading time-scales as the 
viscosity drops sharply during the ring's evolution. 

For a disk in the self-gravitating regime, the viscosity 
is dominated by the translational and gravitational com- 
ponents, and can be expressed as v cx S 2 . In a rough 
description we assume that the disk is a slab of material 
of width L and centred on R = Rq. In that case, the sur- 
face density can be expressed as £ = Mo/ (2ttRqL), where 
Mq is the disk's initial mass. Using L 2 /v cx t we can then 
write L 4 cx t. To check the validity of this relation, the 
disk's width evolution against time is plotted in Fig. \E\ 
We obtain a perfect proportionality relation. 
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Fig. 8. Disk width evolution against time with variable vis- 
cosity. The diamonds are data points from simulations. Note 
how AR 4 remains proportional to time, so that the disk's width 
scales like t 1 ^ 4 . 



3.3.1. Disk mass 

The total mass of the disk is given by: 

M disk = / £ (R) 2ttR(1R 

JO 

On a discrete grid, this equation becomes 

N-l 

Af disk = J2 ^ ( R i+i - R l 



(16) 



(17) 



where N is the number of bins in the grid. 

In Fig. |H]is plotted the evolution of the disk's mass over 
5 billion years, with the CV and NCV models. The NCV 
case can be divided into 3 successive steps. First the disk's 
mass remains constant as no material has yet reached the 
domain's limits. Second, at 2 x 10 7 years the disk reaches 
the domain's outer limit and its mass starts to decrease. 
As seen in the previous section, the disk spreads with a lot 
of material being stacked close to the stiff outer edge. As 
a result, a lot of mass is rapidly lost by the disk when it 
reaches the domain's outer limit. Third, at ~ 10 9 years, 
the disk reaches the domain's inner limit and the mass loss 
is increased. 



1 cr 



' c 



' c 



10 14 



- o 1 



NCV model 
CV mode 



10= 



10° 



10 10° 
Time (years) 



10* 



Fig. 9. Disk mass as a function of time, for a disk with 
variable (solid line) and constant (dashed line) viscosity. The 
disk with constant viscosity is emptied in ~ 10 9 years, while 
a third of the initial mass remains at 5 Gyrs in the disk with 
variable viscosity. 



For a disk with constant viscosity, we directly get L 2 cx 
t. Inverting this relation we get that the width of a disk 
with constant viscosity increases as t 1 / 2 , while with vari- 
able viscosity it increases as t x l 4 . The spreading time-scale 
of a disk with variable viscosity is then much smaller than 
that of a disk with constant viscosity. 

3.3. Mass evolution 

Once the disk reaches the domain's limits (either inner 
or outer), it empties out. We describe below the emptying 
time-scale and show that it is strongly dependant on the 
disk's initial mass and self-gravitating regime. 



While the disk is completely emptied in a few 100 mil- 
lion years in the CV case, the disk mass is still ~ 1 x 10 19 kg 
after 5 billion years in the NCV case. It seems then possi- 
ble for planetary rings to survive viscously over the age of 
the Solar System in the dynamical environment of Saturn 
when considering a physically realistic viscosity model. 

3.3.2. Initial mass influence 

One can argue that the survival of ~ 10 19 kg of mate- 
rial at the end of the disk's evolution is only valid for the 
set of initial parameters considered here, in particular the 
initial mass of 3.75 x 10 19 kg. With a lower initial mass, 
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there would be less material to evacuate from the disk. 
But it would also decrease the surface density and the vis- 
cosity, and thus would slow down the loss of material. The 
opposite behaviour would occur with a higher initial mass. 
Then it raises the question: how does the disk's final mass 
depend on the initial mass ? 

Today's mass of Saturn's rings have been evaluated to 
be si milar to Mimas' mass, or probably more ( Esposito 
et al.. 1 1983T ) . Even if meteoritic flux is expected to bring 
a signi ficant amount of material to the disk f Cuzzi and 
Estrada, 1998h . it seems natural to assume that the rings 
initial mass was larger than their present mass. Anyway, to 
address this question we performed a series of simulations 
starting with initial masses ranging from 1 to 10 Mimas 
masses. Results are plotted in Fig. [TUJ 
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Fig. 10. Disk mass as a function of time, for different disk 
initial masses. With a higher mass, the loss-of-mass rate is 
increased initially but drops down over time, so that all disk's 
mass are comparable after 5 Gyrs of evolution. 



Self- gravitating regime. Initially the disk is fully self-gra- 
vitating. As showed in section 3.2.5 we can then write that 
L 2 oc (M / {2ttR L)) 2 t. Finally we get L 4 cx M§t. 

We can use the above law to estimate the mass lost 
by the disk. Assuming that the disk fills the entire area 
available for its spreading, we get L = constant and the 
further "widening" of the rings is equivalent to a loss of 
mass, because all material passing through the domain's 
boundaries is removed from the disk. Using the above 
expression we can express a time-scale for the emptying of 
the disk t cmpty oc 1/(M§). 

To compare this law with our simulations, the time 
needed for the disk's mass to be half the initial mass, 
against l/Mg, is plotted in Fig. [IT] for initial masses of 1, 
2, 5 and 10 masses of Mimas. We obtain a relation of pro- 
portionality that confirms the above analytical study: in 
the self-gravity regime, the disk is emptied on a time-scale 
inversely proportional to the square of its initial mass. Is 
this relation still valid for a disk in the non-self-gravitating 
regime ? 




The disk that starts with 10 Mimas masses (i.e. 3.75 x 
10 20 kg ) loses 3 x 10 20 kg in the first 100 million years: 
that's 80 % of its initial mass ! So the initial mass of 
the disk does modify the spreading speed and the amount 
of material ejected: adding more mass to the initial disk 
enhances the viscosity and a large amount of material is 
then ejected from the disk very rapidly. In particular, it is 
interesting to note that in order to form rings with a mass 
close to Mimas' mass (i.e. ~3x 10 19 kg), they should not 
be older than 1 billion years old, starting from 10 Mimas 
masses. 

On a larger time-scale, the loss rate drops. Interest- 
ingly, the emptying rate of all disks slows down signifi- 
cantly during the ring's evolution. After 4 billion years, 
no further significant evolution is seen, and at 5 billion 
years, all disk masses are ~ 1.5 x 10 19 kg. 

3.3.3. Analytical insight 

We investigate whether asymptotic laws can be used 
to describe the evolution of the mass lost by the disk, in 
the self-gravitating and non-self-gravitating regime. 



Fig. 11. Timescale for a 2-fold decrease of the disk mass. 
The diamonds are data points from simulations. Note that it 
is proportional to 1/Mo, where Mo is the initial disk mass. 



Non-self-gravitating regime. In the NSG regime, the grav- 
itational viscosity is and the collisional and translational 
viscosities have comparable magnitudes (see Fig. [JJ. Their 
dependence on £ is v co \\ = r 2 i7/35] with /3 = 3/ (4r p p p ), and 
u tTms = 7/311/(1 + /3 2 £ 2 ) with 7 = 0.46(77(20). The total 
viscosity is then 



^coll ~r ^trans — 



AY, + C*£ 3 
1 + /3 2 Z 2 ' 



(18) 



where A 



(3 (r 2 Q + 7) and C = r 2 f2/3 3 . Using the same 
simple approximation 5] = Mq>/ (2-kRqL) and using L 2 jv = 
Kt (K is a proportionality constant) we get 



L 5 + aL 
z{b 2 L 2 + 1) 



Kt, 



(19) 
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y/A/C{2nR /M ) 2 and 



where a = (PM /(2ttR )) 2 , b 
c = C{M /(2nR )) 3 . 

For standard parameters of Saturn's rings (L ps 10 4 km, 
R w 10 5 km, r p = 1 m, p p = 900 kg.m~ 3 ) we get 
L 5 w 10 35 , ai « 10 21 , 6 2 L 2 w 10 2 . We can then sim- 
plify Eq. (pD as: 

Finally the emptying time-scale verifies i cmp ty oc I/Mq ^> 
1 /Mq . As a result the emptying time-scale is much larger 
in the NSG regime than in the SG regime. In the following 
section we investigate the transitions between these two 
regimes, in terms of disk's mass. 

3.3.4- Mass of a fully non self-gravitating disk 

As the disk spreads, the surface density drops down and 
larger parts of the disk (near the edges) become marginally 
self-gravitating. We should then see the emptying time- 
scale drop significantly over time. Fig. [TU] shows indeed an 
initial rapid emptying of the disk, which then considerably 
slows down. All disk's mass even seem to converge toward 
a limiting mass. 

We have shown in the previous paragraph that the 
emptying time-scale is much increased when the disk is 
in the NSG regime. In order to find the transition be- 
tween the SG and NSG regimes, we look for a condi- 
tion on the disk's mass for a disk to be self-gravitating 
or not. We investigate this condition in the case where the 
disk fills the entire area available for its spreading (that is 
from the planet's radius to the Roche limit). We then use 
the Toomre's Q parameter to describe the self-gravitating 
regime of the disk. 

This parameter is a measure of the gravitational stabil- 
ity of the disk. As mentioned earlier we consider here that 
the transition between the SG and NSG regimes occurs for 
Q = 2. The disk is then fully non-self-gravitating if : 



VR,Q 



3.36GS 



> 2. 



(21) 



Using il — y/ GM/R 3 and c r = J Gm p /r p we find that 
the disk is completely NSG if the surface density verifies 



Vi?,£<0.15J^ir 3 / 2 . 



(22) 



Then using (|16p and for the specific case of Saturn, we 
find that the disk is completely in the non-self-gravitating 
regime if its total mass verifies 



M, 



disk 



< M NSG ~ 1-032 x 10 19 kg, 



(23) 



In the simulations presented in this section, the disk's mass 
at t=5 Gyrs is ~ 1.34 — 1.46 x 10 19 kg for initial masses 
of 1 to 10 Mimas' masses. Theses masses are larger than 
Mnsg which should indicate that they are still in the self- 
gravitating regime (at least partially). This is indeed the 
case as seen in Fig. [5] 



In conclusion, it seems that whatever the initial mass, 
the disk undergoes an initial rapid evolution when it is 
self-gravitating. This evolution continuously slows down 
as larger parts of the disk become marginally self-gravita- 
ting. Eventually, the disk will becomes entirely NSG, and 
will undergo a ver y slow evolution. Th e wakes observed 
in Saturn's A ring ( Ferrari et al. . 2009f ) indicate that the 
ring is still self-gravitating, at least in some regions. 

3.3.5. Mass fluxes through boundaries 

In this section, we consider again the previous standard 
model with constant and variable viscosities, starting with 
1 mass of Mimas. We investigate the mass lost by the disk 
through the inner and outer domain's limits, and whether 
theses two quantities are comparable or not. This can 
give clues on the mass falling onto the planet or crossing 
the Roche zone, the latt er being for ins tance available for 



accretion into satellites (jCharnoz et all In press) 



At each time step, the mass fluxes at the limits of the 
domain are tracked. Results are plotted in Fig. rT2J The 
dashed line, representing the mass flow for the CV model, 
is truncated at ~ 1 billion years of evolution, because at 
that point the disk was almost emptied out of its material 
(see Fig. [9]), which caused the numerical code to stop. 
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Fig. 12. Mass passing through domain boundaries (cumula- 
tive) with variable (solid line) and constant (dashed line) vis- 
cosities. The data for the disk with constant viscosity is trun- 
cated at t=1.5 Gyrs, as the disk was empty at that time (see 
Fig. [9j). Whatever the viscosity model, the disk loses much 
more mass through the Roche limit (bottom graph) than by 
infall onto the planet (top graph). 
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Fig. 13. Disk surface density at different times and for different particle radii: O.Of m (black), f m (light blue), 2.5 m (green), 
5 m (dark blue), fO m (red), (a) After f Myrs of evolution, (b) After 100 Myrs of evolution, (c) After 1 Gyrs of evolution, (d) 
After 5 Gyrs of evolution. The disks with larger particles enter more rapidly the non-self-gravitating regime. A disk still partially 
self-gravitating has a specific shape with an inward peak and an outer plateau corresponding to the non-self-gravitating area. 



NCV and CV models behave very differently. For the 
CV disk, both limits of the domain are reached at about 
the same time, ~ 10 6 years. For the NCV disk, both 
limits are reached much later and not simultaneously: the 
outer limit is reached after 20 million years, while the inner 
limit is reached only after 800 million years. This is due 
to the much larger viscosity far from the planet, as the 
translational component scales like i? 5 / 2 in the self-gravi- 
tating regime. 

Another feature, not visible in the log-log format, is 
the very sudden release of material through the outer limit 
(due to the stiff outer edge, see section 3.2.4) compared to 
a smoother and more progressive release through the inner 
edge (due to the smooth inner edge). 

4. Influence of ring parameters 

In order to validate our results of section 3, we now 
perform additional simulations, varying the disk parame- 
ters. In the following, only the variable viscosity model 
(see section 2.2) is considered. 

4-.1- Influence of ring particles radius 

The viscosity model, and our 1-dimensional approach 
requires the fixing of the value of the particle size. In 
this section, we study the influence of this parameter. We 
perform several simulations of the "standard model" , with 
single particle sizes ranging from 0.01 to 10 metres, and 



investigate the impact on the evolution of the surface den- 
sity and the disk's mass, and on the mass fluxes through 
the boundaries. 

Jj.A.1. Influence on surface density evolution 

For a NSG disk, the viscosity increases with particle 
radius (see section 2.2.2.). For a SG disk, the translational 
viscosity does not depend on the particle radius. However, 
the particle radius increases Toomre's Q parameter. As a 
result the disk is more likely to become NSG with bigger 
particles in regions with lower surface densities, i.e. close 
to the edges. Thus, different spreading regimes may then 
be expected whether the disk remains SG during spreading 
or not. In Fig. Q2]the surface mass density of the disk at 
different times, and for several particle sizes, is plotted. 

At 1 million years (Fig. \Wk). the surface density is not 
significantly affected by particle size, except the disk with 
10 m-radius particles (red curve). The ramps close to the 
disk extremities are due to transitions from the SG regime 
to the NSG regime (see section 3.2.4). 

At 10 8 years (Fig. [T3b) the disk with 10 m-radius 
particles has evolved very differently from the other disks. 
This disk has a parabolic shape whereas the disks with 
smaller particles have a surface density peak inward and 
lower densities in the outer regions. This latter shape is 
specific of a disk that is still partially self-gravitating (in 
the peak region). 

At 10 9 years (Fig. [13h). the disk with 5 m-radius 
particles is also completely in the NSG regime and its 
shape is very similar to the disk with 10 m-radius par- 
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tides. The disk with 2.5 m-radius particles is in the SG 
regime only in the region located between the discontinu- 
ity at ~ 75, 000 km and the beginning of the outer ramp 
at ~ 95, 000 km. The surface density of each disk, nor- 
malized by the r p , is plotted on Fig. [T?J along with the 
transition between the SG and NSG regimes (Q = 2), also 
normalized to r p (black dashed line). 




Distance to Saturn (1000 km) 

Fig. 14. Surface density at t — 10 9 years normalized by 
the ring particle radius for different particle radii: 1 m (light 
blue), 2.5 m (green), 5 m(dark blue), 10 m (red). The case 
r p — 0.01 m has been remove for scaling purpose. The black 
dashed line is the transition at Q — 2 between the self-gravita- 
ting to the non-self-gravitating regime normalized by the ring 
particle's radius. The disks with lm and 2.5 m-radius particles 
are still partially in the self-gravitating regime (in the peak), 
while the disks with 5 m and 10 m-radius particles are entirely 
in the non-self-gravitating regime. 

The dark blue and red curves, corresponding to the 
disks with 5 m and 10 m-radius particles, are below the 
dashed line, indicating that they are entirely no longer in 
the SG regime. The light blue and green curves, corre- 
sponding to the disks with 1 m and 2.5 m-radius particles, 
are still partially in the SG regime, in the regions where 
the surface density is above the black dashed line. At 
5 x 10 9 years (Fig. IT3tO. only the disks with 0.01 m-radius 
and lm-radius particles are still partially self-gravitating. 

Two spreading regimes clearly appear. In particular, 
the inner peak of the surface density is characteristic of 
a disk that is still partially self-gravitating, whereas a 
parabolic shape is characteristic of a fully not-self-gravita- 
ting disk. 

Increasing the particles size slows down the initial outer 
spreading, but when all disks are completely NSG, the 
disks with the larger particles spread faster, as viscosities 
in the NSG regime are an increasing function of r p . As a 
result, Fig. [T3H can be interpreted as different evolution- 
ary steps in the disk's life. 

4- 1.2. Influence on disk mass evolution 

The final disk mass for all particle radii is indicated 
in Table [TJ The final mass of the disk is not much af- 



Table 1 

Influence of particle radius on the evolution of the disk's 
mass evolution. Initial mass is 3.75 x 10 19 kg 



Particle 


Final 


Mass through 


Mass through 


radius 


disk mass 


inner edge 


outer edge 


(m) 


(l0 19 kg) 


(10 19 kg) 


(10 19 kg) 


0.01 


1.45 


0.46 


1.84 


1 


1.34 


0.61 


1.80 


2.5 


1.63 


0.61 


1.51 


5 


1.46 


0.69 


1.60 


10 


0.93 


0.93 


1.89 



fected by the change of particles size, except the disk with 
10 m-radius particles because of the strong enhancement 
of the translational and collisional viscosities in the NSG 
regime. The disk with 2.5 m-radius particles is intermedi- 
ate between a disk whose viscosity remains high because 
it remains in the SG regime (for r p < 1 m), and one that 
rapidly enters the NSG regime, but whose viscosity is sup- 
ported by large particles (for r p > 5 m). 

The inner mass flux increases with the particle radius 
because of the increase of the collisional viscosity. The 
outer mass flux increases when the disk enters rapidly the 
NSG regime (r p > 5 m), because of the increase of the 
translational viscosity with r p in the NSG regime. Here 
again the intermediate state of the disk with 2.5 m-radius 
particles can be pointed out. 

4-2. Influence of the initial surface density profile 

In this section we study the influence of different initial 
surface-density profiles. We modify one parameter at a 
time, and analyse its influence on the viscous spreading 
over 5 billion years, and on the evolution of the disk's 
mass. 

4.2.1. Influence of initial mean radius 

First we change the initial position of the disk by shift- 
ing the surface density maximum. While our standard 
model is defined with a mean radius of 110,000 km, we 
study two other cases: 90,000 km and 130,000 km, which 
more or less covers the total area of today's rings. 

The surface density at t = and at 5 billion years is 
plotted in Fig. [15J The final disk masses are listed in 
Tabled 

At 5 billion years, the disks starting at 90,000 km and 
110,000 km have reached a similar surface-density profile, 
and comparable disk mass. The mass fluxes through the 
disk's boundaries have opposite behaviours: while the disk 
starting at 90,000 km loses its mass preferentially through 
its inner boundary, the disk starting at 110,000 km loses 
quite the same amount of mass but through the outer 
boundary (Table 0). 

The disk that starts at 130,000 km reaches the outer 
boundary very rapidly, in about 10 5 years, and loses a 
larger amount of mass than the other two. This is due 
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Fig. 15. Influence of initial position on surface density evo- 
lution. Top: initial profiles. Bottom: at t = 5 Gyrs. When 
the initial ring is farther from the planet it loses more mass, as 
it is then closer to the Roche limit (see Fig. \12} . 



to the fact that when the disk is SG, the translational 
viscosity is greatly increased at high radius, and also to 
the fact that the disk starts closer to the Roche limit. As 
a result its final surface density is globally lower than that 
of the other two disks. This disk never reaches the planet 
within 5 billion years. 

Table 2 

Influence of initial disk position on the evolution of the 
disk's mass. Initial mass is 3.75 x 10 19 kg 



Initial 


Final 


Mass through 


Mass through 


position 


disk mass 


inner edge 


outer edge 


(km) 


(10 19 kg) 


(10 19 kg) 


(10 19 kg) 


90,000 


1,39 


1,58 


0,78 


110,000 


1.34 


0.61 


1.79 


130,000 


0.92 


4 x 10" 4 


2.83 



Eventually, all disks still have a mass close to 10 19 kg, 
much like in the standard model of section 3. The main 
effect of changing the disk's initial position is to modify 
the material fluxes through the boundaries, and increasing 
the initial radius increases the mass loss through the outer 



4-2.2. Influence of initial width 

In this section we change the initial width of the disk by 
adjusting a, the standard deviation of the Gaussian. We 
compare the results of our standard model to simulations 
with standard deviations 0.5a, 2a and 5a. 

The modification of the initial width of the disk has 
no significant impact on the evolution over 5 billion years. 
The most narrow disk, that starts with the highest surface 
density, spreads more rapidly than the others and catches 
them up in a few ten thousand years of evolution. 

4-2.3. Influence of initial mass repartition 

Here we study a scenario where the initial material is 
distributed in two ringlets instead of one. We assume a 
Gaussian profile for both of them and use initial positions 
symmetrical with respect to R = 110,000 km. We choose 
to put the exact same mass in both ringlets, so their initial 
maximal surface density is different (because the surface 
of a bin increases with R). The two sub-rings are distant 
of 10,000 km (peak to peak). The first million years of the 
evolution of the disk is plotted in Fig. [TCI 
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Fig. 16. Surface density evolution over the first million year 
for an initial ring composed of two ringlets (plain line). Both 
ringlets contain the same mass, but the outer one has lower sur- 
face densities because the ring's surface (~ 2nRdR) increases 
with distance. Both ringlets spread similarly (dashed line). 
The two sub-rings merge in ~ 10 6 years (dotted line), after 
which the spreading continues as in the case with a single ini- 
tial ring. 

The two sub-rings merge in a few hundred thousand 
years, and then the evolution goes on very similarly to the 
case with only one initial ring. The final mass of the disk 
after 5 billion years of evolution is 1.34 x 10 19 kg, very 
close to what we obtained in section 3. The mass passing 
through the boundaries is unaffected. The initial shape 
of the rings seems then to have very little effect on their 
long-term evolution. 
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5. Summary and discussion 

5.1. Viscous age of the rings 

We have investigated the effect of a non-constant re- 
alistic viscosity model on the global viscous spreading of 
dense planetary rings over 5 billion years. Our main re- 
sult is that using a physically realistic viscosity model in- 
troduces important changes in the disk's spreading over 
time. We identified two distinctively different spreading 
regimes, whether the disk is self-gravitating or not. When 
the disk is self-gravitating, it develops a surface density 
peak inward and an outward tail with lower surface den- 
sity and which is marginally gravitationally stable (Q ~ 2), 
whereas a totally non self-gravitating disk has a parabolic 
profile with a central maximum. The disk's spreading is 
only significantly affected through modification of the par- 
ticle radius, because bigger particles increase the Toomre's 
Q parameter. 

Contrary to a continuous spreading of the disk with 
constant viscosity, resulting in an emptying in ~ 500 Myrs 
(depending on the constant viscosity value), the spreading 
rate of the non-constant viscosity (NCV) disk significantly 
decreases over time, allowing for a survival of an impor- 
tant part of the disk over 5 billion years. We have shown 
in particular that with variable viscosity in the self-gra- 
vitating regime, the disk's width increases like i 1 / 4 , much 
slower than in a constant viscosity case where the width 
increases like t 1 !" 1 . 

The final state of the disk seems somewhat indepen- 
dent of various initial parameters such as initial width, 
position, and surface density profile. Moreover, it appears 
that whatever the initial mass of the disk, the disk al- 
ways undergo a rapid initial emptying, which progressively 
slows down. We showed that this evolution can be related 
to the disk's self-gravitating regime, with emptying time- 
scales proportional to 1/Mq in the self-gravitating regime, 
and 1/Mq in the non-self-gravitating regime. The disk's 
mass drops continuously until it reaches the mass of a disk 
filling the area from the planet's radius up to the Roche 
limit, and with a surface density so that every part of the 
disk is marginally gravitationally stable (Q ~ 2), yielding 
£ < 0.15-\/-Mmp/rpi?~ 3 / 2 . We computed that for the spe- 
cific case of Saturn this mass is ~ 10 19 kg, in good agree- 
ment with our simulation results (Fig. ITOl) . Afterwards, 
the disk evolves very slowly because of the very low vis- 
cosity in the NSG regime, compared to the SG regime, for 
meter-sized particles. 

In all our simulations starting with a disk mass about 
one Mimas' mass or more, we could not achieve a total 
emptying of the rings in 5 billion years. Most of the mass 
is lost through the Roche limit, except when the initial disk 
is located below ~ 100, 000 km, where most of the mass 
is lost via infall onto the planet. The survival of Saturn's 
rings against viscous spreading over the age of the Solar 
System seems then possible, and during its evolution a 
lot of material would have been made available for the 



The rings of Jupiter, Uranus and Neptune are today 
much less dense and massive than Saturn's rings. More- 
over, they are composed mainly of dust, while Saturn's 
rings have typical particle sizes of ~ 1 cm. For instance, 
the normal optical depth of the Jovian ring system is 
~ 10~ 5 , while the optica l depth of Saturn's A -ring is ~ 0.1 
( Ockert-Bell et all P9991 : iFerrari etHL 12009? ) . In our sim- 
ulations, we showed that under viscous spreading an ini- 
tially massive ring system remains massive over time be- 
cause of the very slowly evolving asymptotic regime when 
the disk becomes marginally self-gravitating. Under these 
considerations, it seems unlikely that the rings of the other 
giant planets could have been dense enough to be self gra- 
vitating (Q < 2) in their past. 

5.2. Comparison with observation 

It is interesting to compare the current rings of Saturn 
with the surface-density profiles we obtain in our simula- 
tions. We have plotted in Fig. [T7]the ring optical depth 
r from the Voyager photo-polarimeter spectrum (PPS), as 
it may be proportional to the surface mass density in low 
t regions. We have also plotted the disk surface density 
at t = 1.5 X 10 s years for the disk with 5 m-radius par- 
ticles. Apart from the Cassini Division, which is thought 
to have been cre ated by Mimas' 2:1 resonance f Goldre- 
ich and Tremaine, 1978h and the C ring, we note several 
similarities. 
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Fig. 17. Voyager PPS profile (plain line) and disk surface 
density at t=1.5 x 10 s years for r p — 5 m (dashed line). Both 
shapes are similar, with a peak inward and lower densities out- 
ward. 



Firstly, the global shape is in good agreement, with a 
maximum inward and lower surface densities toward the 
outer limits of the rings, with the peak corresponding to 
the B ring and the outer plateau to the A ring. This 
plateau is marginally gravitationally stable (Q ~ 2), which 
is co herent with the seve ral observations of wakes in the A 
ring (jFerrari et al. . 20091 ). The position of the maximum is 
much more inward later in our simulations at ~ 10 9 



accretion of satellites (jCharnoz et all In press) 



years. 

This favours the hypothesis that Saturn's rings are young, 
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Fig. 18. Current rings surface density evolution, (a) Initial profile, (b) At 4 Myrs of evolution, (c) At 600 Myrs of evolution, 
(d) At 1 Gyrs of evolution. The B ring rapidly fills the Cassini Division, and overwhelms the C ring in ~ I Gyrs. As in the 
previous simulations, the disk evolves with a peak inward and a marginally self-gravitating outer plateau (see Fig. |6j). 



but this could also be explained by the imperfections of 
the viscosity model we use (the mass has been transferred 
inward too rapidly), or also to other physical effects we 
have not considered in our study, or to unknown initial 
conditions. A change in the particle radius would also 
modify significantly the optical depth profile. 

We can note also that in our simulations we never re- 
cover anything similar to the C ring. This could be due 
to an important modification of the particle size in this 
region, but this is unlikely as particle size distributions 
show a greater power index for the C-ring than for the 
B-ring, so that there shoul d be less large particle s in the 
C-ring than in the B-ring (jCharnoz et ah . 2009al . Section 
17.2.1.1). Moreover this modification of the ring parti- 
cle size would not explain the very steep transition from 
the B ring to the C ring. This transition may however 
suggest that the C ring is a product of another physical 
process. For instance, meteoritic bombardment is though 
to be capable of proje cting large amounts of ma terial over 
significant distances (jCuzzi and Estrada! . 1 199ST ) . The Co- 
ring could then be ejecta from the bombardment of the 
B-ring. This would in turn agree with the observed simi- 
larities in the spectral sign atures of the B and C- rings as 



observed by Cassini VIMS (jNicholson et all 2008) 



Although today's rings mostly resemble a 100 million- 
years old disk when compared to our simulations, their 
formation so recently in the history of the Solar System 
is yet to be explained. In this work we have considered 
only the viscous torque, so that we can only say that 
today's rings viscously look like a 100 million-years old 



disk. Other important physical effects have to be included 
in order to fully constrain the evolution time-scales of 
the rings. For instance, a satellite located between the 
planet and the initial rings could importantly slow down 
the inward spreading of the rings via resonant interactions, 
which appl y a repulsive torque to the disk ( Meyer- Vernet 
and Sicardy, 1987J). Conversely, outer satellites could also 
substantially reduce the rings outward spreading. These 
interactions would consequently reduce the mass lost by 
the disk through its outer boundary. Considering that 
in our simulations, most of the mass is lost through this 
boundary, satellite interactions appear as a process that 
could significantly increase the survival of the rings. How- 
ever these interactions should increase the mass lost by 
the disk from its inner boundary, because of the repulsive 
torque applied by the satellite on the disk. Inclusion of 
satellite interactions would then be important to quantify 
the impact on the mass lost by the disk. 

5. 3. Viscous spreading of current rings 

While the aim of our work was not to reproduce today's 
rings but to understand the effects of a realistic viscosity 
model over long time-scales, we applied our model to an 
initial set of parameters that roughly reproduces today's 
rings of Saturn: a first slab from 74,000 km to 92,000 km 
with S = 100 kg.m~ 2 (the C ring), a second slab from 
92,000 km to 117,000 km with £ = 1000 kg.m" 2 (the B 
ring), a gap from 117,000 km to 122,000 km (the Cassini 
Division) and a third slab from 122,000 km to 136,000 km 
with E = 400 kg.m -2 (the A ring). We choose a particle 
radius of lm. Results are plotted in Fig. [THJ 
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In this configuration, the total rings' mass is 2.2 x 
10 19 kg and only the A and B rings are self-gravitating, 
the C ring being not dense enough. The average initial vis- 
cosities are 10 -4 m 2 .s _1 for the C ring, and 10~ 2 m 2 .s _1 
for the B and A ring s, which is in good ag reement with 



observational results (jTiscareno et all 120071 ). 



The Cassini Division is filled in about 4 million years, 
because of a rapid outward spreading of the B ring (Fig. 
[TBbh This emphasizes the role of the 2:1 Mimas resonance 
in the existence of the Cassini Division (ISchwart A Il98ll : 
Esposito et al.l . ll987T ). The inner spreading of the B ring is 



much slower. The C ring, because of its very low viscosity, 
evolves very little. It is completely "eaten" by the B ring 
in approximately 600 million years (Fig. 118b ). 

It is interesting to note that the shape of the surface 
density at 10 9 years is very similar to what we obtain 
with our standard model (Fig. IT3b-d. light-blue curve, and 
I18t , even though the initial surface density profiles where 
very different in the two simulations. The final state of the 
rings seems then very independent of the initial conditions. 

After 1 billion years of evolution, the mass of the disk 
is still 1.6 x 10 19 kg. It seems then, as we have seen in 
section 3, that Saturn's rings could very well be viscously 
old, but they are also not likely to disappear in the next 
billion years. 

5.4- Perspectives 

An important improvement to this work would be to 
include a particle size distribution in the viscosity model. 
Indeed, we have seen in section 4 that only the ring-particle 
radius significantly affects the spreading of the rings as it 
impacts the positions of the transitions between the self- 
gravitating and non-self-gravitating regimes, for a given 
surface density. This means that if particle segregation 
occurs in the rings, this could lead to modification of the 
viscous behaviour of some regions of the rings. 

While viscous spreading is a key physical phenomenon 
for the evolution of Saturn's rings, it will be affected by 
other effects such as resonant torques with satellite. Many 
structures of the rings are known to be the result of res- 
onan t interactions with the satellit es of Saturn ( Esposito 



et al.. ll987f ). The repulsive torque exerted by outer satel- 



lites could also sign ificantly alter the spreading of the disk 



( Salmon et all . I2009T ) 



These structures, for instance density waves, locally 
modify the density and thus the viscosity. How material 
is transported through these waves could gives us clues on 
the suggested "recycling" of the rings, which could be re- 
sponsi ble for the apparen t low m eteoritic pollution of the 
rings ( Cuzzi and Estrada . 1 1998h . It would also be impor- 
tant to constrain the long term influence of these resonant 
interactions, and see if they increase the possible lifetime 
of Saturn's rings or if they contribute to their destruction. 
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